#!/usr/bin/env python
# coding: utf-8

# ## 农田耕种面积
# version 2021.03.31

# In[ ]:


import cv2
import numpy as np
import json
from math import radians, cos, sin, asin, sqrt, floor, ceil, pi


# In[ ]:


def haversine(lon1, lat1, lon2, lat2): # 经度1，纬度1，经度2，纬度2 （十进制度数）
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    """
    # 将十进制度数转化为弧度
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])

    # haversine公式
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a))
    r = 6371.009 # 地球平均半径，单位为公里
    return c * r * 1000

def distance(point1, point2):
    """计算两个GPS坐标的球面距离，单位：米"""
    return haversine(point1['lng'], point1['lat'], point2['lng'], point2['lat'])



def cv2_to_gps(base_gps: "(latitude, longitude)", xy: (float, float), shift=0, scale=1):
    """计算CV2坐标系中的点的GPS坐标
    @param base_gps CV2的坐标原点的GPS坐标
    @param xy       CV2中待转换的坐标点
    @param shift    CV2坐标相对原点的偏移
    @param scale    CV2坐标的缩放比例

    @return (latitude, longitude)
    """
    x, y = (xy[0] - shift)/scale, (xy[1] - shift)/scale
    base_lat, base_lng = radians(base_gps[0]), radians(base_gps[1])
    r_lat = base_gps[0] + 180 * y / (pi * 6371009)
    r_lng = base_gps[1] + 180 * x / (pi * 6371009 * cos(base_lat))
    return (r_lat, r_lng)



def get_gps_ranges(traces):
    """获取GPS轨迹的经纬度区间范围
    @return 经度区间, 维度区间
    """
    lng_range, lat_range = [180, 0], [90, 0]
    for x in traces:
        lng_range[0] = min(lng_range[0], x['lng'])
        lng_range[1] = max(lng_range[1], x['lng'])
        lat_range[0] = min(lat_range[0], x['lat'])
        lat_range[1] = max(lat_range[1], x['lat'])
    return lng_range, lat_range


def gps_to_cv2(traces, till_width, scale, memory_limit):
    """
    @param traces  GPS数组，每个元素为dict包含'lng'和'lat'
    @param memory_limit 内存限制，单位MB
    """
    lng_range, lat_range = get_gps_ranges(traces)
    # 坐标原点的GPS坐标
    src_pt = {"lng": lng_range[0], "lat": lat_range[0]}
    # 右下角的GPS坐标
    right_bottom_pt = {"lng": lng_range[1], "lat": lat_range[1]}
    # 矩形的高和宽，单位：米
    right_bottom_y = distance({"lng": lng_range[1], "lat": lat_range[0]}, right_bottom_pt)
    right_bottom_x = distance({"lng": lng_range[0], "lat": lat_range[1]}, right_bottom_pt)

    # 自适应scale以便不超过内存限制
    scale = min(scale, sqrt(memory_limit * 1024 * 1024 / (right_bottom_x * right_bottom_y)))
    scale = ceil(till_width * scale) / till_width
    adj_till_width = int(till_width * scale)

    if adj_till_width % 2 != 0:
        adj_till_width += 1
        scale = adj_till_width / till_width

    origin_shift = int(adj_till_width * 0.5 + 1)

    max_point = ((origin_shift + scale * right_bottom_x),  (origin_shift +scale * right_bottom_y))

    results = []
    for g in traces:
        xref_pt = {"lng": lng_range[0], "lat": g["lat"]}
        x = distance(g, xref_pt) * scale + origin_shift
        yref_pt = {"lng": g["lng"], "lat": lat_range[0]}
        y = distance(g, yref_pt) * scale + origin_shift
        results.append((x, y))
    return {
        "points" : results,
        "origin_gps": src_pt,
        "max_point": max_point,
        "scale": scale,
        "origin_shift": origin_shift,
        "adj_till_width": adj_till_width
    }



# In[ ]:


# import psutil

# def get_increased_mem(org_mem):
#     pc_mem =psutil.virtual_memory()
#     return (pc_mem.used - org_mem)/1024/1024


def get_traces_boundaries(traces, till_width=2.5, tolerance=0, scale=5, memory_limit=4000,
                          want_few_points=True, debug=False):
    """获取轨迹覆盖的面积及边框的GPS坐标
    @param traces     GPS坐标数组，每个元素包含'lng'和'lat'
    @param till_width 耕种宽度，单位：米
    @param tolerance  容忍面积参数，当空洞面积小于该参数时忽略空洞，单位：平方米，默认为0
    @param scale      建议缩放比例，实际使用时会考虑内存限制
    @param memory_limit 算法使用的内存限制，单位：MB
    @param debug      是否返回 img(耕种覆盖图在cv2坐标系下的图像) 和 traces(cv2坐标下的边界点坐标列表，
                        列表长度与 boundaries相同)，默认为False

    @return 一个dict对象，包含
            area:       耕种总面积，除去未耕种的空洞，单位：平方米
            boundaries: 边界数组，每个元素包括
                points:  gps坐标数组
                type:    'filled'或'blank'
                subarea: 该边界覆盖的面积，单位平方米
            img:    耕种覆盖图在cv2坐标系下的图像，仅debug=True时返回
            traces: cv2坐标下的边界点坐标列表，列表长度与 boundaries相同，仅debug=True时返回
    """
    if type(traces) is str:
        if traces.find('"')<0: traces = traces.replace("lat", '"lat"').replace("lng", '"lng"')
        traces = json.loads(traces)
#     adj_till_width = int(till_width * scale)
#     origin_shift = int(adj_till_width+2)

    # scale * till_width 至少要大于4
    scale = min(max(scale, 4.0 / till_width), 100 / till_width)

    traces_in_cv2 = gps_to_cv2(traces, till_width=till_width, scale=scale, memory_limit=memory_limit)

    scale, adj_till_width, origin_shift = traces_in_cv2['scale'], traces_in_cv2['adj_till_width'], traces_in_cv2['origin_shift']

    base_gps = traces_in_cv2['origin_gps']
    base_gps = (base_gps['lat'], base_gps['lng'])
    int_pts = np.array(traces_in_cv2['points'], np.int)

    # cv2的坐标顺序是(height, width),返回的坐标点是(x, y),对应关系是x -> width, y -> height
    max_point = traces_in_cv2['max_point']
    img_size = (round(max_point[1] + origin_shift), round(max_point[0] + origin_shift))

    img_cv = np.zeros(img_size, dtype=np.uint8)

    img = cv2.polylines(img_cv, [int_pts], isClosed=False, color=(255,255,255), thickness=adj_till_width)

    # 返回点数由多到少依次为：CHAIN_APPROX_NONE，CHAIN_APPROX_SIMPLE, CHAIN_APPROX_TC89_L1，CHAIN_APPROX_TC89_KCOS
    contours, hierachy = cv2.findContours(img, mode=cv2.RETR_CCOMP, method=cv2.CHAIN_APPROX_TC89_KCOS)

    # 分析区域的类型，计算面积
    boundaries, area = [], 0
    traces = []

    def convert_cv2_to_gps(contour):
        result = []
        for addr in contour:
            gps = cv2_to_gps(base_gps, addr[0], shift=origin_shift, scale=scale)
            result.append({"lat": gps[0], "lng": gps[1]})
        return result

    def refine_points(pts):
        changed = True
        while len(pts) > 3 and changed:
            i, changed = 1, False
            while i + 1 < len(pts):
                if distance(pts[i], pts[i+1]) <= till_width:
                    changed = True
                    pts.pop(i)
                else:
                    i += 1
        return pts

    for ci in range(len(contours)):
        subarea = cv2.contourArea(contours[ci]) / (scale*scale)
        if subarea > tolerance:
            traces.append(contours[ci])
            points = convert_cv2_to_gps(contours[ci])
            if want_few_points:
                points = refine_points(points)
            boundary = {
                "points": points,
                "type": "filled" if hierachy[0, ci, 3] == -1 else "blank",
                "subarea": subarea
            }
            boundaries.append(boundary)
            if boundary['type'] == 'filled':
                area += subarea
            else:
                area -= subarea

    if debug:
        return {
            "area": area,
            "boundaries": boundaries,
            "traces": traces,
            "img": img
        }
    else:
        return {"area": area, "boundaries": boundaries}


# In[ ]:


import sys
import os

if __name__ == '__main__':
    if len(sys.argv) < 4:
        print("""Require at least three arguments:
              1st=json file name,
              2nd=till width,
              3rd=tolerance to ignore hollow,
              4th(optional)=memory limit in MB""")
        sys.exit(1)
    a = sys.argv[1]
    if os.path.isfile(a):
        with open(a) as f:
            traces = f.readline()
    else:
        traces = a
    till_width = float(sys.argv[2])
    tolerance = float(sys.argv[3])
    memory_limit = 4000
    if len(sys.argv) > 4:
        memory_limit = int(sys.argv[4])
    result = get_traces_boundaries(traces = traces, till_width=till_width, tolerance=tolerance,
                                   scale=1, memory_limit=memory_limit, want_few_points=True, debug=False)

    print(result)

